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Type-dependent stochastic Ising model describing the 
dynamics of a non-symmetric feedback modulef 

Manuel Gonzalez-Navarrete^ 


Abstract 

We study an alternative approach to model the dynamical behaviors of biological feedback 
loop, that is, a type-dependent spin system, this class of stochastic models was introduced by 
Fernandez et. al m, and are useful since take account to inherent variability of gene expression. 
We analyze a non-symmetric feedback module being an extension for the repressilator, the hrst 
synthetic biological oscillator, invented by Elowitz and Leibler [7]. We consider a mean-field 
dynamics for a type-dependent Ising model, and then study the empirical-magnetization vector 
representing concentration of molecules. We apply a convergence result from stochastic jump 
processes to deterministic trajectories and present a bifurcation analysis for the associated 
dynamical system. We show that non-symmetric module under study can exhibit very rich 
behaviours, including the empirical oscillations described by repressilator. 


1 Introduction 


Elowitz and Leibler j7j addressed the design and constrnction of a synthetic network providing 
the basic fnnctionality of generating oscillations. Snch fnnctionality is essential to organize time- 
modnlated biological fnnctions like, for instance, the reqnired periodic adjnstment of an organism’s 
physiology to the circadian rhythm [121129]. Their idea was to constrnct a negative feedback loop 
composed by three transcriptional repressors that are not part of any natnral biological clock. The 
resnlting oscillating modnle was called repressilator (see Fignre 1(a) for illnstration). 

We remark, as done by Elowitz and Leibler [7|, that the repressilator displayed noisy behavior, 
this fact mnst be related with stochastic flnctnations of its components. Indeed, Elowitz et al. [S] 
verihed that gene expression is inherently variable, or noisy, dne to random flnctnations in individnal 
cells (see also, Shahrezaei and Swain 1311 and Swain et al. 1331). 

Progress in the nnderstanding of this modnle, as well as others natnrally occnrring networks, 
and their associated control mechanisms demand the development of mathematical models that 
manage good balance between simplicity and nsefnlness. In the literatnre there exist several different 


* Final version, as will appear in Mathematical Biosciences and Engineering 

Tnstituto de Matematica e Estati'stica, Universidade de Sao Paulo. Rua do Matao, 1010, 05508-090, Sao Paulo, 
Brazil, e-mail: manuelg@ime.usp.br 




approaches to this aim. For instance, in the case of the repressilator we conld hnd works snch as 
Chen and Aihara jl], Chen et al. [5] and Wang et ah |37|. 

This kind of interaction modnles are widely known as feedback loops, which carry ont specihc 
fnnctions in a cell, decomposing realistic cellnlar control networks mm- The stndy of feedbacks 
loops have received extensive attention in recent years [31121 EH]. Particnlarly, the development of 
general frameworks can be fonnd in Along |2|, Tyson and Novak [331, and especially in Wang et al. 
|38j . In general, we are interested in nnderstanding for the dynamical behavionr of the concentration 
of molecnles involved in the interactions. 

In this work, we propose an alternative approach, which consists in the application of a class of 
interacting particle systems (see m- The ideas are taken from the type-dependent stochastic spin 
system, proposed by Fernandez et al. IE]. Essentially, the main featnre of this approach is taking 
acconnt the inherent random variability in gene expression. The potential applications in general 
biological signaling networks are discnssed in Fernandez et al. [T3] . 

As an illnstration of onr approach, we propose the stndy of a particnlar cyclic-interaction loop, 
that we call non-symmetric clock module, which is a simple extension of repressilator. Thns, to 
address this kind of qnalitative approach, we propose the application of a mean-held type-dependent 
Ising model dynamics. Althongh the Ising model was initially stndied to nnderstand the physical 
phenomenon of ferromagnetism iig. nowadays this model represents a nsefnl tool in different areas 
snch as image processing, nenral networks or earthqnake dynamics [a 1301135], among others. 

In onr modeling scheme, the dynamics of the type-dependent Ising models is projected onto 
associated continnons-time jnmp processes, called density-profile processes, which are random walks 
mapping the macroscopic evolntion of the particle systems. 

Fernandez et al. [T3| showed that for arbitrary bnt hxed time intervals, in the limit of a very 
large nnmber of particles (thermodynamics limit), the evolntion of these jnmp processes converges 
to time dependent fnnctions satisfying a correspondent deterministic dynamical system. We will 
inclnde a simpler and straightforward proof of the convergence from the stochastic trajectories of 
the density-prohle to deterministic paths rnled by non-linear differential eqnations. Onr techniqne 
is based on the work of Ethier and Knrtz m, who characterized a class of Markov jnmp processes 
called density dependent population processes (see also Knrtz [T9]). 

Therefore, we nse the convergence resnlt to stndy the dynamical behavionrs of onr non-symmetric 
clock modnle. Particnlarly, we will analyse the inflnence of parameters in the behavionrs of onr model. 
The characterization of the evolntion of the associated dynamical systems, that is a bifnrcation 
analysis, completes onr work. 

We remark that the approach introdnced by Fernandez et al. na allows ns to analyse general 
classes of feedback loops. For instance, we can ennmerate the class of feedback loops stndied in IZH, 
which are simple enongh to be analysed in a theoretical framework, bnt admiting very rich dynamical 
behavionrs. Moreover, all the loops have been fonnd in transcriptional regnlatory networks (Leite 
and Wang [21]). 

As we will see in section H the approach proposes a non-reversible stochastic microscopic dy¬ 
namics based on interacting particle systems and statistical mechanics ideas. We nsed a mean-held 
version as a snitable approximation of other Ising-type interactions. In this sense, Mendonga and 
de Oliveira [23] discnssed properties of the stationary measnre of a TDSIM with near neighbonrs 
interaction in the stndy of repressilator. Farther investigations conld focns in other feedback loops, 
exploiting applications and to better nnderstanding of the type-dependent stochastic dynamics. 
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The paper is organized in the following way. In section we introdnce the non-symmetric clock 
modnle, which will be nsed as an illnstration. Section [^inclndes the dehnition of type-dependent spin 
models, and particnlarly the dynamics of the mean-held type-dependent stochastic Ising model. We 
constrnct the associated density-prohle process and prove convergence to deterministic trajectories 
in section]^ Sectioninclndes a bifnrcation analysis of the dynamical system associated to onr non- 
symmetric clock modnle. Finally, section contains some conclnsions and fnrther investigations. 


2 Feedback loops: the example of a non-symmetric clock 
module 

In this section we introdnce a specihc feedback modnle, which will be used to explain the step by 
step of our modeling scheme. However, we remark that the stochastic approach in Section can be 
applied for dynamical studies of general feedback loops. 

We use the notion of network motifs to design the repressilator and our non-symmetric clock 
module. A feedback loop motif is a simple representation of a transcriptional regulation network m 
l36] . That is a cycle in a directed graph whose vertices (also called nodes) can represent concentrations 
of proteins or genes (see Figure [^. In this sense, the edges can represent either positive or negative 
interactions. In other words, a positive (resp. negative) edge implies that the component in tail 
vertex activates (resp. represses) the transcription of the gene of the element in head vertex. 

We shall consider a simple feedback loop that we call non-symmetric clock module, it is based 
on the repressilator proposed by Elowitz and Leibler [7j. The repressilator is a three transcriptional 
repressor systems that are not part of any natural biological clock, and were used to build an oscil¬ 
lating network in Escherichia coli. In other words, this is a negative feedback loop which provides 
the basic functionality of generating oscillations. Which is biologically required in several contexts 
like in cell-cycle and in the setting up of circadian cycle. 


Figure 1(a) shows the representation of the repressilator. We denote by A, B and C, the Lad 
protein, the TetR gene and the Cl gene, respectively. Accordingly, the hrst repressor protein. Lad 
from E. coli, inhibits the transcription of the second repressor gene, tetR from the tetracycline- 
resistance transposon TnlO, whose protein product in turn inhibits the expression of a third gene, 
cl from 1 phage. Einally, Cl inhibits lad expression, completing the cycle. In other words, the rate 
of change of component A density at each time depends only on C density in an inhibitory manner: 
the density of A tends to decrease if the concentration of C is high. A similar dependence holds 
between C and B and between B and A. 

In this work, we propose the study of a non-symmetric cyclic-interaction module, see Figure l(b)[ 
We will refer to A, B and C as being three types of molecules. We focus on this simple feedback 
loop because it is large enough to admit complex dynamical behaviours. The model includes the 
parameter J that measures the strength of the interaction and a parameter 5 G [0,1], which allows 
us to distribute the interactions between each pair of the three components {A,B,C). Note that 
the interactions between neighbour components could be non-symmetric, but the global interaction 
holds the invariance. 

The particular cases when 6 is equal 0 or 1, are similar to the repressilator above dehned. Al¬ 
though, it is important to stress that for positive or negative values of J, we have inhibition or 
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KigUrG li Representation of cyclic feedback modules: (a) Repressilator: a simple model of transcriptional regulation representing the 
loop of feedback inhibition among three components indicated by A, B and C. The blunt arrows indicate inhibition, (b) Motif of our 
non-symmetric clock module (including the rates of interaction). The parameter <5 E [0,1], and for positive or negative values of J, we have 
inhibition or activation cycles, respectively. 


activation cycles, respectively. We are interested in a characterization of dynamical behavionr for 
this modnle, as a function of parameter J and 5. 

Therefore, we describe a microscopic stochastic approach to derive an associated dynamical sys¬ 
tems, which only seeks to incorporate the essential qualitative information about the biochemical 
interactions. This approach is based in the work of Fernandez et ah [13], which borrows ideas from 
interacting particle systems ra. in such a way that spins represent the internal states of components 
of the feedback loop. 


3 Modeling setup: microscopic type-dependent dynamics 

In this section we describe the type-dependent stochastic spin systems proposed by Fernandez et 
al. [I3| to study signaling biological networks. We focus our explanation by studying the dynamical 
behaviours of the non-symmetric clock module exposed in previous section. In particular, we will 
define a type-dependent stochastic Ising model (TDSIM). That is, a microscopic model with Ising- 
type interactions and having dynamical evolution defined by the type-dependent dynamics proposed 
in Fernandez et al. [T3] . 


3.1 A family of interacting particle systems. 

The type-dependent stochastic spin models are a family of stochastic spin-flip systems introduced in 
Fernandez et al. [T3|. This family extends the usual dehnition of particle systems 1221 , to allow 
asymmetric dependence of rates on the energy function, the Hamiltonian. As a consequence of this 
asymmetric dependence, a particularity of these models is its non-reversible stochastic dynamics. 

Next, we explain these models. We need to introduce a set of spin types T, of cardinal fc, 
representing genes or proteins (for instance: Lad, TetR and Cl, in repressilator above). Also, a 
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vertex set V of a simple finite graph of order N = |V|, denoting spatial positions available for each 
one of the types. We call the ordered pair {i,n) E A = T x V a site. Moreover, the set of internal 
states for each i E T, will be denoted by Si = {ai,..., Os.}. Hence, the spin system has site-space A 
and conhguration space Ha “ rijer Sf. For a conhguration a G Ha we denote by a{i,n) the value 
of the spin at site (f,n) (see Figure [^. 



Figure 2; Lattice representation of the site-space A for a type-dependent spin system. Each vertical line represents a type i £ T. At 
horizontal lines we see the spatial positions V. The central spin shows the notation cr(i,n) S Si, for a given configuration a S 


The continuous-time evolution of these models is governed by a non-reversible Glauber spin-hip 
stochastic dynamics. In other words, we have a stochastic evolution in continuous time, for which 
only one particle hips at each transition. Moreover, the dynamics is non-reversible with respect to 


the Gibbs measure (dehned in (3.6)). 


The corresponding rates are determined in terms of a function : Ha —t 
the spin system. To dehne the Hamiltonian, we denote 


the Hamiltonian of 


C = {i = {i, a) : i E T, a E , 


(3,1) 


then, the Hamiltonian of these models is determined by a family of interaction matrices Sn,i '■ C xC —>■ 
M, one for each pair of spatial positions n,l E V. Due to the original applications in signaling 
biological networks, these matrices are not assumed to be symmetric. Particularly, we say 

that a); (j, 6)] indicates the strength of the inhuence that a spin at a site (i, n) G A in internal 

state a E Si has upon a spin at (j, 1) E A that is in internal state b E Sj. 

As we illustrated in the repressilator module (see Figure [I(^ , type A components act upon type B 
components, while the reciprocal interaction does not occur. Then, it is reasonable that the quantity 
Jn,z[(ba); may be noticeably different that the one in the opposite direction Jz,n[(j, ^); (bo)]- 

However, we are not assuming asymmetry with respect spatial positions, because the interaction is 
only associated to types, that is 
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Jn.i [(i, a); U, 6)1 = Ji.n |(i, a); (j, 6)], 


(3.2) 


for n, I G V, i, j E T, a G Si and b E Sj. Consequently, as usual in statistical mechanics, the 
Hamiltonian is then dehned by 

HA{a) = -Y^ (3.3) 

(i,n)eA (i,0eA 

for each conhguration a E fl a- 


3.2 Type-dependent stochastic dynamics. 

We dehne the stochastic dynamics proposed by Fernandez et ah [TB]. This dynamics allows only 
single spin-flips, more precisely, single-site internal-state transitions. In other words, we assume that 
at each time only one molecule could be produced or degraded. In statistical mechanics notation, 
we say that we have a Glauber type stochastic dynamics. 

Formally, given a conhguration a E tiA, a site {i,n) and an internal state a E Si, we denote 
the conhguration with 




a, = ihn), 

otherwise. 


(3.4) 


That is, a conhguration for which we hx the value a for the spin at site Thus, the energy 

cost for the transition from to is given by 


Aa-^-b, 


<^) = ffEL,) - HK.n))- 


(3.5) 


Usually, this total change of the energy associated to the hip at site {i, n) from state a to h, is 
used to dehne the rates of transition in particle systems. Hence, generally we have a continuous-time 
stochastic evolution being reversible with respect to the Gibbs measure. 


hA(t^) = 






(3.6) 


for each cr G Ga- 

However, in type-dependent spin models the dehnition of the rates of transition is quite diherent. 
The asymmetry of the interaction dehned above leads naturally to the decomposition of (3.5) in the 
following manner. 




a) = A[IN]f4l 


a 


a^b, 


+ A[0UT]|,„, 




(3.7) 


where. 
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(3.8) 


(i,OeA 

that collects the change in the inflnence of the conhgnration a npon the site {i,n), when internal 
state there changes from a to b. On the other hand, 

A[OUT]^-5(a) = (3.9) 

OiOeA 

collects the change of the inflnence that the site {i, n) has on all other sites when its internal state 
flips from a to b. 

Fnrthermore, as a particnlarity of the type-dependent spin models, each transition rate depends 


only on the energy changes bronght npon the site (3.8). Thns, we denote A“r^^(cr) the rate of a 
transition flipping to which depends only on ( |3.8| ), that is. 


= 4> (AlIN]J;J5((r)) , 


(3.10) 


where <h is a non-increasing M+-valned fnnction satisfying the detailed-balance condition <h(A)e^ = 
<h(—A)e“^. 

Finally, we state a formal dehnition of the spin models proposed by Fernandez et ah 


Definition 3.1. A type-dependent stochastic spin model is the continuous-time process (crt)t>o on 
defined by a spin model with a type-dependent interaction, given by the Hamiltonian in (3.3) and 


a dynamics with rates of the form (3.10). 


The dynamics of these spin systems yields non-reversibility with respect to the Gibbs measnre 


(3.6). Then, the type-dependent dynamics has mathematical interest for itself. It is natnral to stndy 


the stationary measnre for particnlar cases of the interactions matrices That is, for instance 

local interactions like nearest neighbors (n.n), as usnal in many particles systems. In Mendonga and 
de Oliveira [23], was proposed an Ising-type n.n interaction to analyse, by Monte Carlo simnlations, 
some properties of the stationary measnre for the repressilator. 

In this work, as done in Fernandez et ah [13], we consider a mean-held interaction for an Ising 
model. The mean-held interaction is a natnral approximation of local interactions and given onr 
interest in concentrations of molecnle, as global observable, this is a snitable choice. Of conrse, the 
convergence resnlts to be explained in Section hold whenever we assnme mean-held dynamics. 

Definition 3.2. A type-dependent stochastic spin model is mean-field if the Hamiltonian parameters 
are of the form 


Jn,i[(ba); U^b)] = 


|V| ’ 


(3.11) 


where {cxijjijeC; a real matrix. 
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We conclude this section explaining the particular modeling setup that we will use to study the 
non-symmetric clock module. 


3.3 Mean-field TDSIM 


We follow the ideas of Fernandez et al. [13] to model our non-symmetric clock module through a 
type-dependent stochastic Ising model (TDSIM), which has all internal state spaces Si = {— 
for all i eT. 

We think the spin types as points {A,B,C} over the circle (see Figure l(b)| ) and, for each 
i E T = {A,B,C}, we denote h{i) the neighbour of i in the clockwise direction. And a{i) to be 
the neighbour in the anti-clockwise direction. Borrowing statistical mechanical nomenclature, we say 


that for J < 0, the activation interactions in the Figure 1(b) are ferromagnetic. Moreover, for J > 0, 
the inhibition cycle is said to be an anti-ferromagnetic model. 


Next, observe that the set (3.1) will be particularly dehned by 


C = {(X +1); (A, -1); (B, +1); (B, -1); (C, +1); (C, -1)} . 


(3.12) 


In addition, the set of vertices V will represent a reservoir of capacity N (actually, we have three 
reservoirs, one for each type of molecule in T). In this way, for a conhguration a G { — 1, +1}^, the 
spin value a{i,n) = -|-1 means that there is a molecule of type i E T at spatial position n E V. 
Otherwise, if a{i,n) = —1, it means the absence of molecules of type i G T in the corresponding 
reservoir at spatial position n E V. 


Acordingly, the mean-held TDSIM must be dehned by rates functions as in (3.10), where $(A) = 
3“^ (for details, see [I3])- In addition, the real matrix {aijjij^c* ns in (3.11) is given by 




(3.13) 


—6Jb, if j = h{i) and a = -|-1, 

— (1 — 5)Jh, if i = a{i) and a = -|-1, 

Hi, if j = i, 

0, otherwise, 

6 E [0,1], and Kj G M is an external held (or chemical 

1, will not inhuence others sites. 


for {i,a),{j,b) E C*, and where J E 
potential). It must be clear that a site {i,n) G A in internal state 
because at {i,n) there is no molecule. 

Notice that, as we said above, for 5 = 0 (as well as 5 = 1) and J > 0, the totally asymmetric 


case, we obtain an Ising model for the repressilator represented in Figure l(a 


Finally, the transition rates for TDSIM are dehned for each a E Da by 


(+2 “M ■ -hf + Ki ) , and 


\,n) 


(a) = exp (-2 -jS . at - ^1^ . hf + k , ) , 


(3.14) 


where af = \{l E V : a{a{i),l) = +1}|, and hf = \{l E V : a{h{i),l) = -|-1}|, for i E {A,B,C}. For 
each type i, the parameter Ki must be interpreted as the own rate of production of molecules i eT. 
















In general Ki is a linear function of parameter J. Because we are interested in the relation between 
two-body interaction (J) with individual interaction [k). 

We summarize our modeling setup in the following statement: 


Definition 3.3. The microscopic evolution of the non-symmetric clock module is described by a 
mean-field type-dependent stochastic Ising model with continuous-time Glauber dynamics {o't)t>o, 
whose rates of transition are given by (3.14). 


In the next section we will focus in the macroscopic evolution of the non-symmetric clock mod¬ 
ule. Thus, we shall define the density-profile processes, which are a family of random walks with 
continuous time evolution. Moreover, we will prove an almost sure convergence from these stochastic 
trajectories to deterministic paths governed by non-linear differential equations. 


4 Macroscopic evolution: the density-profile process and its 
convergence to deterministic dynamics 

In this section we focus our attention in the concentrations of biochemical components into our 
feedback loop defined in Section Hence, we will mostly be interested on the vector of empirical 
magnetization of the mean-field TDSIM. Finally, we will study the thermodynamic limit, that is, the 
behaviours when |V| = iV —)■ cx). 

We project the mean-field TDSIM onto a jump Markov process in called density-profile 
process, that represents the densities of the three biochemical components in our feedback loop. 
Thus, in statistical mechanics notation, the stochastic spin system and the density-profile process 
provide, respectively, the microscopic and macroscopic views of the same model. 


Definition 4.1. The density-profile process associated to the non-symmetric clock module is a continuous¬ 
time jump process {X^{t))t>o € = {[0,1] H {k/N,k G for |V| = V > 1. That is defined 

for each t > 0 by 

X^{t)-.= {X^{tfiX^{tfiX^{t)fi (4.1) 

where 


I{n e V : ai(i,n) =+1}| 
- 


(4.2) 


i G {A,B,C}. Moreover, the jumps of this process are of the form \/N, where 1 belong to J = 

{±( 1 , 0 , 0 );±( 0 , 1 , 0 );±( 0 , 0 , 1 )}. 


We describe the evolution of that processes in the following way: for each position X^{t) = a; G 
V]sf, the jumps of the form 


X —)■ X -|- 


1 

iv’ 


(4.3) 
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where 1 G JT”, are defined by a collection (3\{x) of fnnctions, /3i : —)■ [0, cxd], 1 G J. Of course, we 

require that x G and l3i{x) > 0, imply x + \/N G Vj^. It must be clear that each (3\{x) will be 


dehned as a function of the transition rates in (3.14). 


Therefore, the jump in (4.3) occurs with rate 


d 


Np,{x) = -P \X^ {t + s) = X + - 


where 


/3i(x) = 


ds 




N 


(4.4) 


s=0 


if 1 = -e* 
if 1 = e,;. 


(4.5) 


where Xi = {t) as dehned in (4.2), and is the unitary vector in the direction of i G {A, B, C}. 

Note that each variable Xi represents the density of elements of type i, or we say as well, the 
concentration of the component i G {A, 5,0} in our non-symmetric clock module. Thus, Nxi will 
represent the number of molecules of the type i, and in some abstract sense, iV(l — Xi) denotes the 
available number of possible molecules into the reservoir of size N. Although, in the Ising system, 
this last quantity is just the number of sites of type i G T, with spin-value equal to —1. 

The next step is the characterization of these jump processes in the thermodynamics limit, that 
is N ^ oo. Particularly, Fernandez et al. [13], showed that in the limit of a very large number of 
spins, these density-prohle processes converge to time dependent functions satisfying a correspondent 
deterministic dynamical system. They used a pathwise approach, which strongly exploits large 
deviation theory [26|, and coupling of random variables [33]. Their method provides explicit bounds 
for the distance between the stochastic and deterministic trajectories. 

However, the qualitative study of the behaviors of our feedback loop, does not require this kind 
of accurate results obtained in [13]. Therefore, we follow the works of Kurtz [19] and Ethier and 
Kurtz m, to obtain a simpler and straightforward proof of the convergence from the paths of the 
density-profile processes to deterministic trajectories, these latter ruled by non-linear differential 
equations. 

Now, we aim to state our first theorem and prove it. The main tool that we will use is the 
characterization of a family of jump processes given by Kurtz [19]. We stress that our result can 
be easily extended to more general models, that include applications to a very extensive class of 
feedback loops. 

Initially, we define a jump Markov process X^ on for which we will use some properties given 
by Kurtz [TU]. Particularly, the jumps of this process are of the form 1 G (see Dehnition 4.1), and 
their intensities will be given by /3f(/c) = N(3\{k/N), k G Z^, where P\{x) as dehned in (4.5). Thus, 
let X^(0) being non-random, the evolution of X^ at time f > 0, can be written as 




(4,6) 


where Zi{t) := |{s < — X^{s—) = 1}|, that is, the number of jumps of size 1 until time t. 

Then, by some results from Chapter 7 of [13], we have: 
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(4.7) 


x'^it)=x^(o )+^ ly, U/y: . 

for each t > 0, where the hi(-u), 1 G represent independent Poisson processes with corresponding 
intensities u. 

Now, setting 


F(x) = ^lA(x), 


(4,8) 


and = N Thus, we have 


X''(t) = X«(0) + ^y,yj‘l3,{X'^(s))dsyj‘F(X'^(s))ds, (4.9) 

for each f > 0, where Y\{u) = Y\{u) — u, is the Poisson process centred at its expectation. Thus, our 
first result is the almost sure convergence of stochastic trajectories (4.1) to associated deterministic 
paths. 


Theorem 4.2. Consider the density-profile process X^{t), as defined in (4.1)-(4.2), with intensities 
given by (4.4), (4.5) and (3.14), which satisfies (4.9). Suppose limAr_^oo X^ (0) = xq , and X satisfying 


X{t)=xo+ / F{X{s))ds, t>0. 


(4.10) 


Then for every t > 0, 


lim sup |X^(s) — X(s)| = 0 a.s. (4.11) 

s<t 


The proof of this result is essentially repeat the arguments of Theorem 2.1 from Chapter 11 of 
HU, but we include the proof for completeness. 

Proof. First of all, denote := sup 3 ,g 25 ^ fifix). Thus, it is easy to see that 


|l|/?i < oo, 


(4.12) 


and by the Lipschitz continuity of the rates of the jumps (4.5), there exists a fixed M > 0 such that 


\F{x) - F(y)\ < M\x - y\, x.yeVf,. 


(4.13) 


Observe that by (4.9) and (4.10), we have for each f > 0, 
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\X^{t)-X{t)\ < |X^(0)-xo| + 




+ 


[F{X^{s)) - FiX{s))]ds 


(4.14) 


Now, let denote 


£Nit) ■= sup 

U<t 


r p,{x'^(s))ds 

j \ '^ 0 


(4.15) 


Therefore, using (4.13), (4.15) and the Gronwall inequality, we can rewrite (4.14) as follows. 


\X^{t)-Xit)\ < \X^(0)-xo\+s^it)+ M\X^is)-X{s)\ds 

Jo (4.16) 

< (|X^(0)-xo|+£iv(^))e^^ 


for all f > 0. Thus, since e^it) is a non-decreasing function. 


sup |X^^(s) - X(s)| < (|X"(0) - Xol + SN{t))e 

S<t 


Mt 


(4.17) 


To obtain (4.11) we need the following lemma. 


Lemma 4.3. For every f > 0, 


lim = 0 a.s. 

N^OO 


(4.18) 


Proof. First, note that we have the following equivalence in distribution 


Yi{Nu) D ^ Yf{u) 


N 


i=l 


N 


(4.19) 


where T'j*(-u) are independent Poisson processes with rate u. Then, the strong law of large number 
implies that 


lim sup|A^ ^yi(iVM)|=0 a.s., v > 0. 

Al—>00 


(4.20) 


Furthermore, by dehnition of /9i above, it follows that 

Yi (Nu^i) 


^Nit) < ^ |1| sup 


U<t 


N 


(4.21) 


12 

















If u <t, a basic property of Poisson processes states that, 

(ivwA) < A {Nm) , 

for each 1 G Then, 


|1| sup 

U<t 


Yy [Nu^y) 


N 


<^-^{Yy{Nt^)+NtPy) 
, , A {Nt^y) , , - 


N 


for all > 1 and each 1 G JT”. Now, note that. 


1 1 i=l 




N 


i=l 


(4.22) 


(4.23) 


(4.24) 


where Y^{u) are again independent Poisson processes with instensity u. Finally, by the law of large 
numbers applied to the independent random variables in bracket at r.h.s of (4.24), and by condition 

diia, 


lin, W|l d‘P‘ft) ^ 

N^oo ^ N ^ 


|i|tA 




iii lim 

^ AT-s-oo ^ 


2=1 


Yim 

N ■ 


(4.25) 


That is, by (|4.23|) we can interchange the limit and summation for the expression at r.h.s in 

□ 


(4.21). Therefore, using (4.20), the Lemma follows. 


Finally, from inequality (4.17), by condition lim 7 v->.oo X^(0) = xq, and Lemma 4.3, then for every 
t > 0, 


lim sup |X'^(s) — X(s)| =0 a.s. 

N^oo 


(4.26) 

□ 


Next section is devoted to our results about modeling a specihc biological feedback loop, that is, 
our non-symmetric clock module. Clearly, the application of the Theorem 4^ allows us to analyse the 
qualitative behaviour of the concentrations of the molecules involved in the interactions. Therefore, 
we are able to include the role of stochasticity in the gene expression, but we also simplify the 
qualitative study through a bifurcation analysis for the associated dynamical system. 
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5 The associated dynamical system: formulation and bifur¬ 
cation analysis 

In this section we include a qualitative study of the dynamical behaviours of the non-symmetric clock 
module. We use the ideas of previous sections to analyse a dynamical system related to stochastic 
evolution of the concentration of molecules involved in our feedback loop. 

First, we briefly summarize the ideas in previous sections. Initially, we stated a mean-held TDSIM 


{o't)t>o on IIa and dehned its dynamics by the rates of transition (3.14). Thus, the concentrations 
of each component in our non-symmetric loop (types A, B and C) were characterized by a density- 
prohle process, dehned in (4.1). In the Ising case, the only independent variables are the densities of 


activated types that we denoted i G {A, B, C}. The rate of a jump in density-prohle process 

was dehned in (4.4) and (4.5), based on the rate of corresponding transition in TDSIM. After that, 
in Section]^ we also studied the thermodynamic limit of the particle system {N —)■ cx)), proving that 
the density-prohle process converges to deterministic trajectories governed by non-linear diherential 


equations (see Theorem 4.2) 


Finally, in this section we show that depending on the parameter values, the magnetization 
random path can either converges to a unique stable hxed point (if — 1 < J < 2), converges to one 
of a pair of stable hxed points (for J < —1), or asymptotically evolves close to a deterministic orbit 
in (when J > 2). 


Therefore, we follow Theorem 4.2 to study the dynamics of our non-symmetric clock module by 
analysing the associated dynamical system. However, we remark that the behaviours of the dynamical 
system are deterministic approximations of the stochastic evolutions of the jump processes in our 
approach. 

The system of diherential equations associated to the cycle-interaction module in Figure l(b)| 
and with rates given by (3.14) will be expressed by 


X, 


— (1 _ X '(^ (1 


(5.1) 


for i E {A^B,C}^ J G M, 5 G [0,1], and Ki G M. As usual, we will consider Ki = J/2, i G {A, H,^}, 
that it, proportional to the molecules interaction parameter. The particular value is to reduce the 
number of parameters and to obtain suitable hxed points in dynamical analysis (steady state of 
concentrations (1/2,1/2,1/2)). 

Our next result is the statement of a bifurcation analysis to guarantee for the non-symmetric 
clock module that: for 0 < J < Jc = 2, the concentration of all three components {A,B,C) remains 
stable, but oscillates with large amplitude as soon as J increases past the threshold Jc = 2. For 
J < 0, we can see the appearance of two stable points when J < —1. These situations are showed in 
Figure 


Proposition 5.1. 

= J/2, then 


Consider the dynamieal system (5.1), with 6 E D = [0,1/2) U (1/2,1], and 


K 


a) For J < 0, there is a bifurcation at Jc = —1.' the fixed point (1/2,1/2,1/2)'^ loses the stability 
and appear two stable points for J < Jc- 
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X 



J 


KigUrG 3; Bifurcation diagram for the non-symmetric clock module. The vertical axis represents the concentration of one of the three 
molecules. At right hand, we see a Hopf bifurcation with respect to the parameter J, that measures the strength of the interactions. Solid 
lines indicate stable points; dotted lines indicate unstable points; black circles indicate maximum and minimum values of stable orbits; 
bifurcation diagram obtained with xpp-aut HQ]. 


b) If J > 0, there is a Hopf (or Andronov-Hopf) bifureation at Jc = 2. 


Proof. First, it is easy to see that the dynamical system has a hxed point at = (1/2,1/2,1/2)^, 
becanse x = F(x°, J, 6) = (0, 0, 0)^, for all J G M and 6 E D. 

From now on, we consider the dynamical system y = F*{y,J,6), where y = (|/i,?/ 2 , 2 / 3 )^ and 
Xi = \ + yi. Then, the fixed point for the new system is 0 = (0, 0, 0)^. 

Following the theory in Perko 123 and Knznetsov pni. we conld gnarantee that near 0 = ( 0 , 0 , 0 )^, 
the dynamical system is locally topologically eqnivalent to its linearisation y = Ay, where A denotes 
the Jacobian matrix evalnated at 0. Then, y^ is stable if all eigenvalnes A of A satisfy Re A < 0, 
and stability is lost when one of the real parts becomes positive. 

Notice that A = dF*[A) is given by the expression 

/ 1 (1-5)J dJ \ 

-2 5J 1 (1-J)J . 

5J 1 J 

The eigenvalnes of this matrix are 

Ai = — J — 1, 

A 2 = J (1 — 26)\/3i + (J — 2 ), and 
As = J{26 - l)V3i + {J - 2). 

Therefore, the fixed point y^ loses the stability at Jc = —1, when Ai crosses the imaginary 
axis throngh the origin. Thns, we have two stable fixed points. This bifnrcation is called fold (or 
pitchfork) bifnrcation. The valnes of these fixed points could be calculated by solving: sinh(2J|/) + 
2y cosh.{2 Jy) = 0, thus we obtain the relation 


(5.2) 


(5.3) 
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(5.4) 



for 7 ^ 0. We notice the symmetry around J-axis, this fact could be seen at left hand of Figure 


On the other hand, the stability is lost when A2 and A3, symmetric around the real axis, cross the 


imaginary axis. This fact occurs when J = 2. The phenomenon includes the appearance of a stable 


orbit, this kind of bifurcation is known as Hopf bifurcation. 


□ 


The oscillations experimentally verihed in the repressilator (Elowitz and Leibler [7j) are explained 
in this dynamical behaviour through the (supercritical) Hopf bifurcation with respect to the inter¬ 


action strength parameter J. 

5.1 Deterministic macroscopic evolution in details 

Let us hnally show a detailed analysis of the dynamical behaviours of the system using linearisation 


tools. In fact, we do not study the dynamical system (5.1), however we analyse the linear system 


that approximates it. 

We will exhibit a qualitative analysis of the concentration of the molecules involved in the in¬ 
teractions. The idea is a better understanding of the influence of parameter 5 in the dynamical 
behaviours. 

Therefore, we must translate and rotate the coordinate system X 1 X 2 X 3 onto a new system 
Z 1 Z 2 Z 3 . Particularly, this new coordinate system has the origin at (1/2,1/2,1/2). Thus, will be 
the diagonal, from ( 0 , 0 , 0 ) to ( 1 , 1 , 1 ); the Z 2 -axis is directed to x = ( 0 , 1 , 1 / 2 ); and the Zi-axis is 
oriented to x = (3/4, 3/4, 0). 

Equivalently, remember that y = (yi, ?/ 2 , Xi = \ + yi, for i = {1,2, 3}. We can say that 

the system Z 1 Z 2 Z 3 is obtained by the following two operations in the system YiY' 2 ^ 3 : 

1. - Rotate the yiy 2 b 3 -system around the Is-axis by a = 45° (anticlockwise direction); 

2. - Rotate the Yiy 2 F 3 -system again about the new rotated y 2 -axis by /3 = 55° (clockwise direction); 

Finally, to obtain the new coordinates of the system, we have ( 2 : 1 , Z 2 , Z 3 ) = {yi,y 2 , ys) ■ R, with 


1 1 1 

Ve V 2 P3 


(5.5) 


R 


2 n 1 
\ Ve ^ V3/ 


Again, by the theory in [201 [22], the behaviour of the dynamical system in Z 1 Z 2 Z 3 , could be 
studied by linearisation, that is, we analyse 
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where A = R ^ ■ dF*{Q) ■ R, because R is an orthonormal matrix. Then, 


/ J-2 V^J{25-1) 


0 \ 


A= -x/3J(2h-l) J-2 

Vo 0 


0 


(2J + 2)/ 


0 


(5.7) 


and expression (5.6) represents the following system: 



(5.8) 


Notice that the behaviour over Zs-axis is independent of the dynamics on ZiZ 2 -plane. Then, we 
could establish the following two elementary conclusions: 

Cl) For any J > —1, the trajectory on Za-axis goes to zero, because £3 is negative for Z 3 > 0, and 
positive for ^3 < 0; 

C2) For any J < 2, the dynamics on ZiZ 2 -plane goes to zero in both directions. In other words, it 
goes to the point (0,0). Of course, the velocity and direction depend on 5 and J. 

Accordingly, given Cl) and C2), we conclude that for —1 < J < 2 the system has a stable point 


at 0 = ( 0 , 0 , 0 )^. 


Now, we study the dynamics for J = —1. In this situation, £3 = 0, for any Z 3 . Then, in adition 
with C2), any initial condition 2 :( 0 ) = (zi(0), 2 : 2 ( 0 ), 2 : 3 ( 0 ))^ goes to (0, 0, 2 : 3 ( 0 ))^, as t —>■ cx). So, 
for J < —1 it is easy to see that, ia < 0 when 2:3 is negative, and becomes positive when Z 3 > 0. 


Namely, if the initial condition is established above the ZiZ 2 -plane, the trajectory goes to (0, 0, 77)^, 


as t —)■ cx), where K = K{J, 6 ) is a positive constant that depends on J and 6 . When the system 


starts at the bottom of the ZiZ 2 -plane, the trajectory goes to (0, 0, —K)'^, as f —)■ 00 . These are the 


two stable points in Proposition 5.1 


Furthermore, considering J = 2, since Cl), we just need to study the dynamics over ZiZ 2 -plane. 
That is. 



(5.9) 






We could change variables to polar coordinates. Let zi = rcosO, Z 2 = rsinO. To derive a 
differential equation for r, we note z'l^- Z2 = r^, so ziZi + 2:2^2 = rr. Substituting for zi and £2 yields 
r = 0. Thus from 9 = (2:1^2 — Z 2 Zi)/r‘^, we find 6 = —2y/3(26 — 1). Therefore, the origin is a center, 
and any initial condition z{0) = (2:1(0), 2:2(0), ^3(0))^ evolves over the ZiZ2-plane on the circle with 
radius r(0) = 2:1(0)^+ 2:2(0)^, the direction and velocity of the cycle depend on S. The Figure]^ shows 
us the behaviour of the cycle for different values of 6 . 




Figure 4; Phase portraits for J = 2 and different values of 5. At left hand, we see the behaviour for <5 < 1/2. The right picture show 
US the cycles for 5 > 1/2. The highest velocities will be reach at 5 = 0 and 5 = 1. 


In the remaining case, for J > 2, we must study the two-dimensional linear system z = Az with 

/ J-2 V3J{26-1)\ 

A={ (5.10) 

\-y/3J{26-l) J-2 J 

thus, changing variables to polar coordinates we obtain the system r = r(J —2), and 6 = —'/3J{26 — 
1). Notice that, r > 0 for all r > 0, and the value of 9 depends on 6 . See Figure]^ 




KigUrG Si Phase portraits for J > 2. The left picture shows us the behaviours for 6 < 1/2. At right side, we see the cycles for 5 > 1/2. 


6 Conclusions and comments 

On one hand, based on a particular feedback loop we have described a theoretical framework to 
model and analyse the non-linear dynamics of gene regulatory networks. In the literature there exist 
several different approaches to this aim. However, in previous sections we showed that our approach. 
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based on the well studied Ising model, can include the inherent variability in gene expression, also this 
microscopic approach could capture detailed characteristics of the system interpreted as a promotion- 
inhibition circuitry. Thus, it may be particularly useful to illuminate how simple feedback loops 
manage to perform their basic functionalities and how these functionalities are further integrated 
into the whole cellular regulatory network. 

On the other hand, the mathematical study of the type-dependent dynamics suggests a long¬ 
standing theme of statistical physics of nonequilibrium systems. That is, the question of the nature 
of their stationary measure, and in particular of the existence of phase transitions. Thus, it will be 
interesting for us to propose local interactions in type-dependent Ising models, instead of mean-field 
interaction, at this point we could consider several lattices (like or triangular lattice) and treat 
these questions. 

Particularly, the study of the TDSIM with local interactions is closely related with works like 
Godreche [H] and Godreche and Bray [15] (see also de Oliveira [25]). They considered a kind of 
asymmetric Ising dynamics, called directed Ising model. In [TH |T3] was studied several dynamics, 
like Glauber or Metropolis. They obtained a large space of parameters defining the rates functions 
allowing irreversible Gibbsian Ising models, whenever the dynamics is not totally asymmetric. 

Furthermore, bifurcations showed in the associated dynamical systems (Section suggest the 
existence of phase transitions in the type-dependent Ising model. Of course, it would also be natural 
to consider metastability issues [26] to characterize the dynamical behaviors of this particle system. In 
this sense, a related work is due to Kotecky and Olivieri la, they studied a discrete-time Metropolis 
dynamics for a two dimensional ferromagnetic asymmetric Ising model, in which the vertical and 
horizontal interaction parameters have different values. 

Finally, we say that another very interesting mathematical question is the characterization of 
loss and recovery of Gibbsianness in stochastic evolutions, considered in recent works, in the area 
of mathematical physics. For instance, van Enter et ah [9] studied some Ising models for which 
this phenomenon occurs. In addition, Kulske and Le Ny [18] initiated a fruitful research direction 
showing that Gibbs-non-Gibbs transitions can also be defined naturally for mean-field models. 
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